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Abstract. - The nucleation of vapor bubbles within a superheated fluid is studied using density 
functional theory. The nudged elastic band technique is used to find the minimum energy path- 
way from the metastable uniform liquid to the stable uniform gas thus emphasizing the analogy 
between the the nucleation problem and that of chemical reactions. The result is both an accurate 
determination of the critical nucleus and an unbiased description of the density profile at various 
points along the path between the free energy extrema. This calculation is compared to two other 
methods: the use of parametrized profiles and constrained minimization of the free energy. The 
results indicate that the recent claim, based on the constraint method, that bubble nucleation and 
growth involves an activated instability is incorrect. 
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Introduction. — The process of nucleation, e.g. of 
solids from supercooled liquids or of liquids from super- 
saturated gases, has become the subject of intense interest 
in recent years driven by applications in physics, chem- 
istry and biology. Practical problems such as the control 
of protein crystallization [1] have motivated experimental, 
theoretical and computational research which has uncov- 
ered new fundamental physics such as the important role 
that intermediate metastable phases may play in nucleat- 
ing new phases [1—5]. In cases where this is important, 
nucleation becomes more complex than the simple picture 
of barrier crossing that underlies classical nucleation the- 
ory (CNT). Indeed, there is some theoretical evidence that 
the nucleation of solids from gases (a simple model of the 
precipitation of solids from solution) may be nonclassical 
even in the case of simple fluids [5] . Recently, it has even 
been claimed that the nucleation of gas bubbles in a su- 
perheated fluid may be nonclassical due to the presence of 
an "activated instability" [6]. 

Density functional theory (DFT) provides a natural ap- 
proach to the description of nucleation. DFT is based 
on the fact [7, 8] that the free energy is a unique func- 
tional of the local density, p(r). In simple systems, dif- 
ferent phases correspond to different local densities: for 
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example, in gases and liquids, the density is a constant 
p(r) = p whereas in a solid it varies over the length scale 
of the lattice. Nucleation can then be viewed as a process 
of moving from one point in density-function space to an- 
other via the most probable path. This is analogous to 
the usual description of a reaction pathway in a chemical 
system. As is usually done in chemical problems, it is nat- 
ural to assume that the most probably path will be the 
minimum free-energy path (MFEP) so that the problem 
reduces to one of finding the MFEP connecting the initial 
and final states given the model free energy functional. 

The local density is a field so that in principle its value 
at every point in space can vary independently giving it 
an infinite number of degrees of freedom. Here, attention 
will be restricted to spherically symmetric profiles which 
will furthermore be discretized on a lattice giving a finite, 
but still large, number of degrees of freedom. The free en- 
ergy is therefore a surface in this high dimensional space 
making its characterization difficult. This is analogous to 
the problem of characterizing reaction paths in chemical 
systems. Fortunately, recent years have seen the develop- 
ment of several sophisticated methods for addressing the 
problem of finding optimal pathways on high-dimensional 
surfaces(i.e., eigenvalue following, string method, elastic 
band methods ...) [9]. One of the primary goals of this 
paper is to show that these methods can be applied to the 
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description of nucleation. 

Liquid-vapor nucleation in simple fluids has, of course, 
been studied previously and might seem too classical to 
warrant further effort. Although previous work has mostly 
focused on nucleation of droplets in a supersaturated gas, 
the issues are the same for droplet and bubble nucleation. 
Early studies [10, 11] focused on the determination of the 
critical cluster since this is accessible (it is an extrcmum 
of the free energy) and determines the free energy barrier 
and, hence, the nucleation rate. In the 1990's Rciss et 
al [12-14] developed an approximate simulation technique 
for studying non-critical clusters and this led to paral- 
lel theoretical developments [15] which have recently been 
extended to the study of gas bubbles in superheated flu- 
ids [6]. (This is referred to as the "constraint" method 
below.) For use in simulations, this method assumes a 
decoupling of clusters from the surrounding bulk. More 
recently, importance-sampling techniques for Monte Carlo 
simulations have been developed which do not depend on 
assumptions like decoupling [16,17]. These developments 
in simulation techniques are a primary motivation for the 
present work which aims to develop a similarly unbiased 
theoretical description of nucleation. 

In this paper, three methods will be used to determine 
the properties of bubble nucleation. The simplest is the 
use of a parametrized profile. It is assumed that the den- 
sity profile of a bubble is sigmoidal and is approximated 
by a generalization of a hyperbolic tangent. These profiles 
have two important parameters: the location of the inter- 
face, i.e. the bubble radius, and the width of the interface. 
The approximate pathway is constructed by minimizing 
the free energy with respect to the width of the interface 
for different values of the radius. The second method is the 
constrained optimization method originated by Talanquer 
and Oxtoby [15] for the study of droplet nucleation and 
adopted by Uline and Corti (UC) [6] for bubble nucleation 
where the free energy is minimized while subject to a con- 
straint that enforces a particular bubble geometry. The 
third method is the nudged elastic band (NEB) [18] which 
is a chain of states method for constructing the MFEP. It 
will be shown that of the three methods, the NEB pro- 
vides the most natural and reliable determination of the 
reaction path and that the instability identified by UC is 
an artifact of their method. 



DFT calculations. 
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In DFT, the grand potential is 



pil\p]=PF\p]-/3n P {r)dr 



<(r)p(r)dr (1) 



where (3 = l/fc^T, p is the chemical potential, (r) is 
an external one-body field and p (r) is the local density. 
The functional (3F[p\ is related to the Helmholtz free en- 
ergy [8, 19]. To be precise, DFT tells us that for fixed T , 
p and <fi(r), the equilibrium density is a stationary point 
of the functional (30 [p] . In this sense, there is little direct 
relevance of this functional to the problem of nucleation 



since nucleation typically occurs at some given fixed field 
(either no field or the gravitational field) so that this func- 
tional can only be used to describe a few simple states: the 
bulk liquid and vapor and, potentially, the critical nucleus. 
All other intermediate states will, in general, not extrem- 
ize the functional. However, it is usually assumed that at 
fixed field, (30 [p] plays the role of a potential and that 
it can therefore be used to deduce interesting quantities 
such as the energy barrier for transition from one state to 
another. 

The calculations reported here will be based on a re- 
cently developed DFT model consisting of a modified Fun- 
damental Measure model for the short-ranged repulsion 
and a Van der Waals tail for the long-ranged attraction 
[20]. This model gives a quantitatively accurate descrip- 
tion of the inhomogencous Lennard- Jones liquid including 
the planar liquid-vapor interface, the liquid near a wall and 
the liquid confined to a slit pore. The bulk thermodynam- 
ics was approximated using the empirical 33-parameter 
equation of state of Johnson, Zollweg and Gubbins [21]. 

Three methods are used to explore the free energy land- 
scape. The first and simplest is the use of a parametrized 
profile. The parametrization used here is a modified hy- 
pertangent, 

. 1 + br 1 - tanh (Air - R)) ,„, 



l + (brf l-tanh(-AR) 
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and where pi (p v ) is the density of the bulk liquid (vapor) 
for the given temperature and chemical potential. This 
function is constructed so that it approaches pi as 1/r at 
large distances while also satisfying dp/dr — > as r — > 
0. In the calculations presented below, the free energy is 
minimized with respect to the width parameter A for fixed 
values of the radius, R. 

The second method used was the constraint method 
of Talanquer and Oxtoby [15] as adapted by Uline and 
Corti [6]. Here, one characterizes a bubble by some prop- 
erty of the density which can be written schematically as 
g([p],T) = where the notation indicates that the left 
hand side is a functional of the local density and that it 
also depends on one or more parameters denoted collec- 
tively as r. Different values of the parameters correspond 
to bubbles having different properties (mass, radius, ...). 
It is assumed that for fixed parameters, the most energet- 
ically accessible structure is that which minimizes the free 
energy subject to the structural constraint. This is im- 
plemented via a Lagrange multiplier by extremizing the 
functional (30, [p] — ag ([p] , T) giving the equations 
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and g([p] ,r) = 0. In this picture, the role of the con- 
straint is to reduce the number of degrees of freedom from 
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an essentially infinite number needed to describe all pos- 
sible density functionals to a few, the collection of param- 
eters r, which then define a path from the uniform liquid 
to the uniform vapor. The question then is what to use 
for the structural constraint. First, note that if there is no 
structural constraint, then the uniform state satisfies the 
usual thermodynamic relation — 4^- = p, where / (p) is 
the bulk free energy per unit volume, and there are two so- 
lutions corresponding to the uniform vapor, p v , and liquid, 
pi, states. Perhaps the most natural choice for the con- 
straint is the number of atoms " missing" from the bubble 
relative to the liquid background, 



9([p},T) = J( Pl -p(r))dr-AN 



(5) 



where AN is the desired deficit. In this case, the only 
parameter is T = AiV . However, substitution into Eq.(4) 
shows that this changes the chemical potential from p to 
p, + a so that the liquid density far from the bubble would 
necessarily be incorrect. This seems to indicate that a 
localized bubble can only be described by a constraint 
that vanishes at large r . UC therefore propose 



guc([p],T)= / 9(A-r) p(r)dr-N 



(G) 



where 6(x) is the step function and there are two param- 
eters, the "radius" A and the total number of atoms in the 
bubble, N [6]. (For the reverse problem of droplet nucle- 
ation, Talanquer and Oxtoby also use this constraint but 
they say they adjust the bulk outside the volume in a way 
which appears to imply the use of Eq.(5) [15].) UC show 
that, for a fixed value of N, solutions only exist for A less 
than some critical value and this is interpreted as indi- 
cating that bubbles of larger size are unstable. I propose 
here another constraint that is based on defining the size 
of a bubble implicitly as the position at which the density 
becomes greater than some value, p», 



<?*(M,r) 



&(p*-p(r))dr-T. 



(7) 



This is physically appealing and is also similar to a cri- 
terion used to define clusters in some computer simula- 
tions [16]. In this case there are in principle two parame- 
ters: r, which is a measure of the volume of the bubble, 
and p* which will be taken to be h [pi + p v ) . 

Finally, the third method used is the nudged elastic 
band [18]. These calculations begin with some number, 
N, of density profiles corresponding to Eq.(2) with values 
of R spaced evenly from R = to some R ma x- The various 
profiles represent points along a path from the initial state 
(the uniform fluid corresponding to R = in Eq.(2)) to 
the final state which here is taken to be again a sigmoidal 
with a fixed, large radius and optimized width. (These 
endpoints remain fixed during the calculation and the re- 
maining images are referred to as "moveable".) Unlike 
the first method, each of these images are not restricted 



to be given by the parametrized form of Eq.(2) but is in- 
stead a discretized density profile, p(ri), that can assume 
any shape. If the sum of the free energies of all of the 
moveable profiles was minimized, the final result would be 
that each profile would tend toward one of the minima - 
cither the uniform liquid or the uniform vapor, depending 
on the starting point. To avoid this, the NEB involves two 
elements. First is the addition of artificial elastic forces be- 
tween neighboring profiles so that the images are forced to 
be more or less uniformly distributed along the MFEP. To 
define the elastic forces, the distance between two profiles 
was defined to be 



\Pl - P2\ 



w(f)(px{r) - p-zir)) dr. 



(8) 



The definition allows for a weight, w(r), since there are at 
least two natural choices one could make in a spherical ge- 
ometry. The first is, of course, w(r) = 1 while the second 



is w(r) 



For spherically symmetric density distribu- 



tions, the first choice has the effect of giving more weight 
to differences in density at large r than at small distances 
while the second weights all radial points equally. Calcu- 
lations were performed using both definitions and it was 
found, as shown below, that the choice had no effect on 
the final results except for affecting the relative spacing 
of the images along the MFEP: the second choice gives 
more images at small r than the first choice but the path 
so defined was virtually identical. 

The second clement of the NEB method is that this com- 
bination of the free energies and artificial clastic forces is 
minimized under the constraint that the clastic forces only 
act along the MFEP and the force arising from the free 
energies act only perpendicular to the MFEP (the "nudg- 
ing"). This method is not guaranteed to give the exact 
MFEP but it has proven a useful heuristic in numerous 
applications. A further modification involves choosing one 
image and forcing it to climb the gradient along the MFEP 
without any spring forces [22] . This gives a very accurate 
and convenient characterization of the saddle point. The 
relaxation of the system of images and harmonic forces 
was performed using the the fast inertial relaxation en- 
gine algorithm [23]. 

Results. — The three methods were compared 
for the case of a Lennard- Jones potential, V(r) = 

4e ((7) 12 - iff) with a cutoff at r = 4cr. The calcula- 
tions reported here were performed at a reduced temper- 
ature ksT/e = 0.8 and the chemical potential was fixed 
at 10% below coexistence corresponding, roughly, to the 
state conditions used by UC [6] . The densities of the co- 
existing liquid and vapor were found to be pier 3 = 0.784 
and p v a 3 = 0.008 respectively. 

Figure 1 shows the excess free energy, calculated using 
all three methods, as a function of the size of the bub- 
ble as quantified by the number deficit (i.e. as defined 
via eq.(5)). The constraint method has been implemented 
using eq.(7). All three methods are in close agreement 
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Fig. 1: (Color online) Excess free energy as a function of the 
atomic deficit as calculated using DFT with the proposed struc- 
tural constraint given in Eq.(7), circles and solid line, using a 
modified hyperbolic tangent profile, Eq.(2), squares and broken 
lines, and the NEB method, diamonds. 



for the height of the free energy barrier and the size of 
the critical cluster. However, away from the critical clus- 
ter, the constraint method gives systematically higher free 
energies than do the other two methods. For the initial 
part of the bubble formation, from the uniform fluid to 
the critical cluster, the NEB method gives the lowest free 
energies and, hence, is closest to the MFEP. For bubbles 
larger than the critical cluster, the NEB energies are very 
close to, but slightly higher than, the those calculated us- 
ing the parametrized profile. This is probably due to the 
fact that the end point of the chain of states used in the 
NEB calculation is not the uniform gas - since the uniform 
gas is a bubble of infinite size, it cannot serve as the end of 
the chain - but, rather, a parametrized profile which thus 
distorts the NEB chain away from the MFEP. 

Figure 2 shows the density profiles calculated using the 
NEB method. The bubble initially forms as a small drop 
in the density very close to the origin. The central density 
quickly becomes that of the coexisting gas but the bub- 
ble is still very small due to the fact that the interface 
is extremely narrow. At this stage, some structure forms 
in the liquid which is similar to that seen in a fluid near 
a hard wall [20]. The bubble grows via the broadening 
of the interface and as the interface broadens, the liquid 
structure vanishes. Only when the interface is sufficiently 
broad, in this case, when the number deficit is about 70 
atoms or so, does a recognizable sigmoidal structure de- 
velop. Further growth consists of a displacement of the 
interface with little change in structure. 

When the free energy is minimized for a spherically sym- 
metric density profile with the UC constraint, the results 
of UC are reproduced: namely, for fixed N, there is a max- 
imum value of the radial parameter, A, beyond which no 
stable profile exists. The physics of this apparent instabil- 
ity will be the subject of the next section. 

Using the proposed constraint, a different picture 
emerges. It is possible to stabilize the density profile for 
all values of the parameter T. For each value of T, the to- 
tal atomic deficit can be calculated and the resulting free 



0.8 


_ ~J 1 ' I 


- 


0.6 


ffvvv 




" 




-If;/ i / ■' / 


N=0.01 




III/ 8 / / i 


N=0.2 






\\i i / 1 / 


— N=l 




0.4 


1 / / / ' / 


N=3 

N=7 


" 




\ \ § / / / 


— N=13 






III / 1 1 / 


o N=13 
N=73 


"" 


0.2 


\\\ / / / / 


— N=145 
N=280 


- 




yyJJ . 


1 



2 4 6 8 

r/o 

Fig. 2: (Color online) Structure of small bubbles along the 
MFEP as determined using the NEB method. Structures up 
to N = 13 were calculated using the metric with w(r) = r -2 
weight whereas structures for N = 13 and larger were per- 
formed including the w(r) = 1. The N — 13 structure is 
shown as calculated both ways (line is without the weight and 
the symbols are with the weight) illustrating the insensitivity 
of the calculation to the weight used in the metric. 



energy curve is shown in Fig. 1 which confirms the CNT 
picture and shows no sign of an instability. However, be- 
cause the constraint has the effect of shifting the chemical 
potential inside the bubble, but not outside, (see eq.(4)), 
the profiles generated are discontinuous, as can be seen in 
Fig. 3. The only continuous profile obtained is for the 
special case of the critical nucleus, for which the Lagrange 
multiplier a is zero. The discontinuity of the profiles is 
responsible, at least in part, for the higher free energies 
calculated for this model as shown in Fig. 1. 

Illustration on a simple model. It is useful 

to consider a simple, but instructive, analytic model for 
which it is possible to see why the constraint method can 
give instabilities. Consider the square-gradient model first 
introduced by van der Waals [24,25], 



I3F{ P ] 



Pf{p(r)) + \r 



dp(r) 

Or 



dr (9) 



where again / (p) is the free energy per unit volume of the 
uniform fluid and where the constant r is called the accom- 
modation factor. This model is not quantitatively very 
accurate, but it is a systematic approximation to more 
realistic models and is still often used in the literature. 
Combined with the structural constraint, this can be used 
to determine the density. The problem is simplified if one 
uses a parametrized density and here I consider a very 
simple picccwise linear approximation, 



p(r) = 



po, r < R 

Po + £2 it Bi (r-R), R<r < R + w 
Poo, R <r 



(10) 



Potentially, the width, w, the bubble density, po , the outer 
density, poo and radius R are determined by minimizing 
the grand potential (with the structural constraint). For 
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Fig. 3: (Color online) Density profiles for several different bub- 
ble sizes of, from left to right, AN = 4, 163, 288 and 866 as 
calculated using the constraint method. 



the outer density, this gives p^ = pi. However, for the 
remaining parameters, this is still too complicated, so I 
consider the simplest case that the width goes to zero. If 
this limit is taken, while holding 7 = t/w constant, then 
the resulting functional using the UC constraint is easily 
shown to be 



/3fl[p}-ag uc ([p},T) (11) 

Aqr 

= yi? 3 {/(Po)-/(p/)-M(Po-p/)} 

+47T 7 i? 2 (pi - p ) 2 

(f A 3 p - N) Q(R - A) 
(f R 3 Pa + f{\ 3 - R 3 ) Pi - N) 9(A - R) 



Notice that in the absence of the constraint, a = 0, this is 
the usual functional assumed in classical nucleation the- 
ory if one takes po = p v . In that case, the free energy 
is a function of a single parameter, the radius R, and it 
is easy to see that (3£l [p] as a function of R shows a min- 
imum at R = (the liquid) a maximum at some finite 
value of R and a monotonic decrease for larger values of R. 
This is just the CNT picture of nucleation as a process of 
thermally activated barrier crossing. Using the proposed 
constraint, the functional to be minimized is the same as 
the UC functional except the Lagrange multiplier term is 
—a (^fR 3 — r), provided that p < p* < p;. 

The difficulties of the UC model become apparent when 
the free energy is minimized with respect to the free pa- 
rameters of the density profile. There are two cases de- 
pending on whether R or A is larger. Minimizing under 
the assumption < R < A gives 
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47 (pi - poY 
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(12) 



— N - (A 3 - R 3 



4tt 



Pi 



Now, consider what happens when A becomes large and TV 
is held fixed. Since p > 0, it follows that R 3 > A 3 - -r^—N 
which means that R ~ A. This in turn implies that po < 



^ 3 R i N ~ j§^N so that the central density tends to zero. 
However, at small po, one has that / (po) = Po hipo — po + 
O (Po) giving 



R = 



47P? 



O(p ) 



pilnpo - f(pi) + O(p ) 



(13) 



and this is clearly negative for sufficiently small po, which 
is to say for sufficiently large A. Since the radius must be 
positive, this is unphysical. The alternative, R > A, leads 
to 
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Po 



27 (pi - PoY 



f(po) - f(Pl)~ P*(P0 

3 N 



Pi) 
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Again, the central density tends to zero as A becomes 
large. In this case, one has that 



lim R 



ZlPl 2 



HPl 



-f( Pl )+ m (3P{ Pl ) 



(15) 



where P(pi) is the pressure of the bulk liquid. It is evident 
that the radius again becomes negative thus violating the 
assumption that R > X. There is therefore no nontrivial 
solution for large A, which is exactly the conclusion of UC. 
However, in this case the underlying free energy is clearly 
well defined for all bubbles (i.e. for all values of p and 
R) so that the lack of solution is a statement about the 
structural constraint and tells us nothing about the free 
energy of bubbles. 

It might seem that the pathology of this model is due to 
the extremely simplified density profile. However, a look 
at the proposed structural constraint shows that this is 
not necessarily the case. In fact, minimization of the free 
energy with respect to the Lagrange multiplier, a, and the 
central density, po, gives T 



^-R 3 and 



f'(Po) = P 



37 (pi - poY 

R 



(16) 



The first relation shows a simple one-to-one relation be- 
tween the parameter T and the bubble radius while the 
second relation serves to fix the central density. For large 
radii, that density tends to the density of the bulk gas. For 
small radii, taking into account that the chemical potential 
for a gas is generally negative and that f'{po) ~ log(po), it 
is clear that the central density tends to zero. There is a 
third relation that results from minimizing with respect to 
R but this serves to fix the Lagrange multiplier and is of 
no particular interest. Finally, note that these conclusions 
are independent of the value of p* provided it lies in the 
range [p v ,p{\. 

This analysis is consistent with the more detailed cal- 
culations of the previous section and shows that the lack 
of existence of a density profile that minimizes the con- 
strained free energy does not necessarily mean that the 
unconstrained free energy is pathological in any way. In 
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the present case, the CNT free energy is a well-behaved 
function of the local density (characterized by the param- 
eters poj Poo and R). 

Conclusions. — The MFEP for bubble formation in 
a superheated liquid was calculated using three different 
methods. The most robust method appears to be the 
NEB. In this approach, no a priori assumptions are made 
about the structure of bubbles and the results were found 
to be robust with respect to the metric used and the choice 
of parameters. The method yields a well-defined picture 
of bubble formation and growth that is broadly consistent 
with that of classical nuclcation theory. 

It is perhaps surprising that the use of parametrized 
profiles appears to be almost as good as the NEB method 
in terms of the description of the critical cluster as well 
as the calculation of the free energy as a function of num- 
ber deficit. For very small bubbles, the assumption of a 
sigmoidal shape is not correct, but the free energy is rela- 
tively insensitive to this error since, in a spherical geom- 
etry, the characteristic r 2 weighting in the integrals min- 
imizes the contribution of the structure near the origin 
while larger bubbles actually are sigmoidal. 

The constraint method gave the worst results. In par- 
ticular, it is not robust: the use of the UC constraint only 
gives nontrivial results for small bubbles whereas the alter- 
native constraint discussed her gives stable profiles of all 
sizes. In both cases, the profiles are discontinuous except 
for the critical cluster thus giving a higher free energy than 
the other methods. The failure of the UC constraint to 
stabilize large bubbles was analyzed using a simple model 
and it was shown that the apparent instability was an ar- 
tifact of the constraint method and tells nothing about the 
underlying free energy surface. 

The conclusion, based on the combination of numeri- 
cal calculations and analysis of a simple model presented 
here, is that bubble nucleation and growth does not in- 
volve an "activated instability" as previously claimed [6] 
but, rather, appears to follow the picture, if not the de- 
tails, of classical nucleation theory rather well. 
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